## R script for Rost & Clarke. 2023. "Correlates of Aid:
## Why do governments fund some humanitarian appeals but not others?"

## Load required R packages
library (fixest)
library (betareg)
library (corrplot)
library (car)
library (sjPlot)
library (countrycode)
library (lmtest)
library (performance)
library (sandwich)
library (dplyr)

## Two datasets are required for these analyses:
## "aid" (for most models) and "aid05missingadded" for Model 4

## Generate coverage99 variable (coverage limited to 1%-99%) in both datasets
## There are a few appeals with zero or 100% or more coverage, however, 
## ...beta regression requires the outcome variable to be between 0 and 1
## There is one appeal with zero coverage: Honduras in 2015
## There are 14 appeals with 100% of higher coverage, 
## ...up to 179% for a flash appeal in Afghanistan in 2021

aid$coverage99 <- ifelse (aid$coverage < 1, 
                          ifelse (aid$coverage == 0, 
                                  .01, 
                                  aid$coverage), 
                          .99)
aid05missingadded$coverage99 <- ifelse (aid05missingadded$coverage < 1, 
                                        ifelse (aid05missingadded$coverage == 0, 
                                                .01, 
                                                aid05missingadded$coverage), 
                                        .99)

## Model 1: Fixed Effects (country and year) OLS 
## ...with 'amount of funding' as outcome variable
fit_fe_base <- feols ( I(log(funded + 1)) ~ I(log(funded.lag1 + 1)) 
                       + I(log(displaced_conflict + 1)) 
                       + I(log(hosted + 1)) 
                       + famine 
                       + I(log(battledeaths + 1)) 
                       + I(log(people_disaster + 1)) 
                       + trade_open 
                       + colony_top10 
                       + unsc_mtgs 
                       + any_pko 
                       + log(news)
                       | year + iso3, 
                       data = aid)

summary (fit_fe_base, vcov = "twoway")
check_model (fit_fe_base)


## Model 2: Fixed Effects (year) Beta Regression 
## ...with 'percent of appeal covered' as outcome variable

fit_beta_fe_year <- betareg (coverage99 ~ coverage.lag1 
                             + I(log(displaced_conflict + 1)) 
                             + I(log(hosted + 1)) 
                             + famine 
                             + I(log(battledeaths + 1)) 
                             + I(log(people_disaster + 1)) 
                             + trade_open 
                             + colony_top10 
                             + unsc_mtgs 
                             + any_pko 
                             + log(news) 
                             + as.factor(year), 
                             data = aid, 
                             link = "logit")

coeftest (fit_beta_fe_year, vcov. = vcovHAC)
summary (fit_beta_fe_year)
check_model (fit_beta_fe_year)

## Correlation matrix for variables used in the two main models (Models 1 and 2)
## (referenced but graph does not appear in the paper)
correlations <- select (aid, 
                        funded, 
                        coverage, 
                        displaced_conflict, 
                        hosted, 
                        famine, 
                        people_disaster, 
                        trade_open, 
                        colony_top10, 
                        unsc_mtgs, 
                        any_pko, 
                        news)

colnames (correlations) <- c ("Absolute funding", 
                              "Percent covered", 
                              "Displaced people", 
                              "Hosted refugees", 
                              "Famine", 
                              "Affected people", 
                              "Trade openness", 
                              "Former colony", 
                              "UNSC meetings", 
                              "Peacekeeping", 
                              "Media articles")

correlations <- as.matrix (correlations)

correlations <- cor (correlations, 
                     use = "pairwise.complete.obs")

corrplot (correlations, 
          method = 'number', 
          tl.col = "dark grey", 
          tl.cex = 0.8, 
          number.cex = 0.8)


## Check for normality of residuals 
## (referenced but graphs do not appear in the paper)

## Model 1
res_fe_base <- residuals (fit_fe_base)
qqPlot (res_fe_base)

## Model 2
res_beta_base <- residuals (fit_beta_base)
qqPlot (res_beta_base)


## Substantive effects

## Model 2: hosted refugees
plot_model (fit_beta_fe_year, 
            type = "pred", 
            terms = "hosted", 
            title = "Model 2: Predicted appeal coverage per number of hosted refugees", 
            axis.title = c("Number of hosted refugees", 
                           "Predicted appeal coverage")) 

## Model 2: trade openness
plot_model (fit_beta_fe_year, 
            type = "pred", 
            terms = "trade_open", 
            title = "Model 2: Predicted appeal coverage by trade openness", 
            axis.title = c("Trade openness", 
                           "Predicted appeal coverage"))

## Model 2: number of news articles
plot_model(fit_beta_fe_year, 
           type = "pred", 
           terms = "news", 
           title = "Model 2: Predicted appeal coverage by number of news articles", 
           axis.title = c("Number of news articles", 
                          "Predicted appeal coverage"))

## Model 10: magnitude of genocide
plot_model(fit_beta_genomag, 
           type = "pred", 
           terms = "geno_mag", 
           title = "Model 10: Predicted appeal coverage for severity of genocide", 
           axis.title = c("Severity of genocide", 
                          "Predicted appeal coverage"))

## Model 11: human rights violations (not shown)
plot_model(fit_beta_pts, type = "pred", 
           terms = "hrights", 
           title = "Model 11: Predicted appeal coverage level of human rights violations", 
           axis.title = c("Human rights violations", 
                          "Predicted appeal coverage"))

## Robustness Tests

## Model 3: Minimalist (fewer variables, more observations, 
## ...including flash appeals and Covid-19 appeals)
fit_beta_small <- betareg (coverage99 ~ I(log(displaced_conflict + 1)) 
                           + famine 
                           + I(log(battledeaths + 1)) 
                           + colony_top10 
                           + unsc_mtgs 
                           + any_pko 
                           + log(news) 
                           + as.factor(year), 
                           data = aid)

coeftest (fit_beta_small, vcov. = vcovHAC)
summary (fit_beta_small)
check_model (fit_beta_small)

## Model 4: Maximalist (many missing values replaced with zeros)
## This is the only model that uses a different dataset
fit_beta_max <- betareg (coverage99 ~ coverage.lag1 
                         + I(log(displaced_conflict + 1)) 
                         + I(log(hosted + 1)) 
                         + famine 
                         + I(log(battledeaths + 1)) 
                         + I(log(people_disaster + 1)) 
                         + trade_open 
                         + colony_top10 
                         + unsc_mtgs 
                         + any_pko 
                         + log(news) 
                         + as.factor(year), 
                         data = aid05missingadded, 
                         link = "logit")

coeftest (fit_beta_max, vcov. = vcovHAC)
summary (fit_beta_max)
check_model (fit_beta_max)

## Model 5: Donor Fatigue / Proctated Crises
fit_beta_fe_protracted <- betareg (coverage99 ~ coverage.lag1 
                                   + years_hrp 
                                   + I(log(displaced_conflict + 1)) 
                                   + I(log(hosted + 1)) 
                                   + famine 
                                   + I(log(battledeaths + 1)) 
                                   + I(log(people_disaster + 1)) 
                                   + trade_open 
                                   + colony_top10 
                                   + unsc_mtgs 
                                   + any_pko 
                                   + log(news) 
                                   + as.factor(year), 
                                   data = aid, 
                                   link = "logit")

coeftest (fit_beta_fe_protracted, vcov. = vcovHAC)
summary (fit_beta_fe_protracted)
check_model (fit_beta_fe_protracted)

## Model 6: Multiple Appeals
fit_beta_multiple <- betareg (coverage99 ~ coverage.lag1 
                              + multiple_appeals 
                              + I(log(displaced_conflict + 1)) 
                              + I(log(hosted + 1)) 
                              + famine 
                              + I(log(battledeaths + 1)) 
                              + I(log(people_disaster + 1)) 
                              + trade_open 
                              + colony_top10 
                              + unsc_mtgs 
                              + any_pko 
                              + log(news) 
                              + as.factor(year), 
                              data = aid, 
                              link = "logit")

coeftest (fit_beta_multiple, vcov. = vcovHAC)
summary (fit_beta_multiple)
check_model (fit_beta_multiple)

## Model 7: Food Security (crisis and higher levels, IPC3+)
fit_beta_ipc3 <- betareg (coverage99 ~ coverage.lag1 
                          + years_hrp 
                          + I(log(displaced_conflict + 1)) 
                          + I(log(hosted + 1)) 
                          + log(ipc3plus) 
                          + I(log(battledeaths + 1)) 
                          + I(log(people_disaster + 1)) 
                          + trade_open 
                          + colony_top10 
                          + unsc_mtgs 
                          + any_pko 
                          + log(news) 
                          + as.factor(year), 
                          data = aid, 
                          link = "logit")

coeftest (fit_beta_ipc3, vcov. = vcovHAC)
summary (fit_beta_ipc3)
check_model (fit_beta_ipc3)

## Model 8: Food Security (emergency and catastrophic/famine level, IPC4 and 5)
fit_beta_ipc45 <- betareg (coverage99 ~ coverage.lag1 
                           + I(log(displaced_conflict + 1)) 
                           + I(log(hosted + 1)) 
                           + I(log(ipc4 + 1)) 
                           + I(log(ipc5 + 1)) 
                           + I(log(battledeaths + 1)) 
                           + I(log(people_disaster + 1)) 
                           + trade_open 
                           + colony_top10 
                           + unsc_mtgs 
                           + any_pko 
                           + log(news) 
                           + as.factor(year), 
                           data = aid, 
                           link = "logit")

coeftest (fit_beta_ipc45, vcov. = vcovHAC)
summary (fit_beta_ipc45)
check_model (fit_beta_ipc45)

## Model 9: Intrastate conflict instead of battle deaths
fit_beta_conflict <- betareg (coverage99 ~ coverage.lag1 
                              + I(log(displaced_conflict + 1)) 
                              + I(log(hosted + 1)) 
                              + famine 
                              + civilconflict 
                              + civilwar 
                              + internationalized 
                              + I(log(people_disaster + 1)) 
                              + trade_open 
                              + colony_top10 
                              + unsc_mtgs 
                              + any_pko 
                              + log(news) 
                              + as.factor(year), 
                              data = aid, 
                              link = "logit")

coeftest (fit_beta_conflict, vcov. = vcovHAC)
summary (fit_beta_conflict)
check_model (fit_beta_conflict)

## Model 10: Magnitude of genocide
fit_beta_genomag <- betareg (coverage99 ~ coverage.lag1 
                             + I(log(displaced_conflict + 1)) 
                             + I(log(hosted + 1)) 
                             + famine 
                             + I(log(battledeaths + 1)) 
                             + I(log(people_disaster + 1)) 
                             + geno_mag 
                             + trade_open 
                             + colony_top10 
                             + unsc_mtgs 
                             + any_pko 
                             + log(news) 
                             + as.factor(year), 
                             data = aid, 
                             link = "logit")

coeftest (fit_beta_genomag, vcov. = vcovHAC)
summary (fit_beta_genomag)
check_model (fit_beta_genomag)

## Model 11: Human rights (Political Terror Scale, PTS)
fit_beta_pts <- betareg (coverage99 ~ coverage.lag1 
                         + I(log(displaced_conflict + 1)) 
                         + I(log(hosted + 1)) 
                         + famine 
                         + I(log(battledeaths + 1)) 
                         + I(log(people_disaster + 1)) 
                         + hrights 
                         + trade_open 
                         + colony_top10 
                         + unsc_mtgs 
                         + any_pko 
                         + log(news) 
                         + as.factor(year), 
                         data = aid, 
                         link = "logit")

coeftest (fit_beta_pts, vcov. = vcovHAC)
summary (fit_beta_pts)
check_model (fit_beta_pts)

## Model 12: Controls (GDP per capita, population size and democracy)
fit_beta_pts <- betareg (coverage99 ~ coverage.lag1 
                         + I(log(displaced_conflict + 1)) 
                         + I(log(hosted + 1)) 
                         + famine 
                         + I(log(battledeaths + 1)) 
                         + I(log(people_disaster + 1)) 
                         + hrights 
                         + trade_open 
                         + colony_top10 
                         + unsc_mtgs 
                         + any_pko 
                         + log(news) 
                         + as.factor(year), 
                         data = aid, 
                         link = "logit")

coeftest (fit_beta_controls, vcov. = vcovHAC)
summary (fit_beta_controls)
check_model (fit_beta_controls)


## Models referred to but not shown in the paper
## Model 1 without peacekeeping
fit_fe_base_nopko <- feols ( I(log(funded + 1)) ~ I(log(funded.lag1 + 1)) 
                             + I(log(displaced_conflict + 1)) 
                             + I(log(hosted + 1)) 
                             + famine 
                             + I(log(battledeaths + 1)) 
                             + I(log(people_disaster + 1)) 
                             + trade_open 
                             + colony_top10 
                             + unsc_mtgs 
                             + log(news) 
                             | year + iso3, 
                             data = aid)

summary (fit_fe_base_nopko, vcov = "twoway")
check_model (fit_fe_base_nopko)

## Model 1 without UNSC meetings
fit_fe_base_nounsc <- feols ( I(log(funded + 1)) ~ I(log(funded.lag1 + 1)) 
                              + I(log(displaced_conflict + 1)) 
                              + I(log(hosted + 1)) 
                              + famine 
                              + I(log(battledeaths + 1)) 
                              + I(log(people_disaster + 1)) 
                              + trade_open 
                              + colony_top10 
                              + any_pko 
                              + log(news) 
                              | year + iso3, 
                              data = aid)

summary (fit_fe_base_nounsc, vcov = "twoway")
check_model (fit_fe_base_nounsc)

## Model 2 without peacekeeping
fit_beta_fe_year_nopko <- betareg (coverage99 ~ coverage.lag1 
                                   + I(log(displaced_conflict + 1)) 
                                   + I(log(hosted + 1)) 
                                   + famine 
                                   + I(log(battledeaths + 1)) 
                                   + I(log(people_disaster + 1)) 
                                   + trade_open 
                                   + colony_top10 
                                   + unsc_mtgs 
                                   + log(news) 
                                   + as.factor(year), 
                                   data = aid, 
                                   link = "logit")

coeftest (fit_beta_fe_year_nopko, vcov. = vcovHAC)
summary (fit_beta_fe_year_nopko)
check_model (fit_beta_fe_year_nopko)

## Model 2 without UNSC meetings
fit_beta_fe_year_nounsc <- betareg (coverage99 ~ coverage.lag1 
                                    + I(log(displaced_conflict + 1)) 
                                    + I(log(hosted + 1)) 
                                    + famine 
                                    + I(log(battledeaths + 1)) 
                                    + I(log(people_disaster + 1)) 
                                    + trade_open 
                                    + colony_top10 
                                    + any_pko 
                                    + log(news) 
                                    + as.factor(year), 
                                    data = aid, 
                                    link = "logit")
coeftest (fit_beta_fe_year_nounsc, vcov. = vcovHAC)
summary (fit_beta_fe_year_nounsc)
check_model (fit_beta_fe_year_nounsc)

## Model 13: Genocide/politicde as a binary variable instead of magnitude
fit_beta_genocide <- betareg (coverage99 ~ coverage.lag1 
                              + I(log(displaced_conflict + 1)) 
                              + I(log(hosted + 1)) 
                              + famine 
                              + I(log(battledeaths + 1)) 
                              + I(log(people_disaster + 1)) 
                              + genopoliticide 
                              + trade_open 
                              + colony_top10 
                              + unsc_mtgs 
                              + any_pko 
                              + log(news) 
                              + as.factor(year), 
                              data = aid, 
                              link = "logit")

coeftest (fit_beta_genocide, vcov. = vcovHAC)
summary (fit_beta_genocide)
check_model (fit_beta_genocide)

## Model 14: only UN peacekeeping instead of any peacekeeping operation
fit_beta_UNpko <- betareg (coverage99 ~ coverage.lag1 
                           + I(log(displaced_conflict + 1)) 
                           + I(log(hosted + 1)) 
                           + famine + I(log(battledeaths + 1)) 
                           + I(log(people_disaster + 1)) 
                           + trade_open 
                           + colony_top10 
                           + unsc_mtgs 
                           + un_pko 
                           + log(news) 
                           + as.factor(year), 
                           data = aid, 
                           link = "logit")

coeftest (fit_beta_UNpko, vcov. = vcovHAC)
summary (fit_beta_UNpko)
check_model (fit_beta_UNpko)

## Model 15: Regions
aid$region <- countrycode (aid$iso3, 
                           origin = "iso3c", 
                           destination = "un.region.name")

aid$subregion <- countrycode (aid$iso3, 
                              origin = "iso3c", 
                              destination = "un.regionsub.name")

fit_beta_subregion <- betareg (coverage99 ~ coverage.lag1 
                               + I(log(displaced_conflict + 1)) 
                               + I(log(hosted + 1)) 
                               + famine 
                               + I(log(battledeaths + 1)) 
                               + I(log(people_disaster + 1)) 
                               + trade_open 
                               + colony_top10 
                               + unsc_mtgs 
                               + any_pko 
                               + log(news) 
                               + as.factor (subregion)
                               + as.factor(year), 
                               data = aid, 
                               link = "logit")

coeftest (fit_beta_subregion, vcov. = vcovHAC)
summary (fit_beta_subregion)
check_model (fit_beta_subregion)

## Model 16: Binary outcome variable

## Create binary variable "wellfunded" = 1 if coverage >= 75%
aid$wellfunded <- ifelse (aid$coverage >= .75, 1, 0)

## The same for the lagged outcome variable
aid$wellfunded.lag1 <- ifelse (aid$coverage.lag1 >= .75, 1, 0)

fit_logit_fe <- glm (wellfunded ~ wellfunded.lag1 
                       + I(log(displaced_conflict + 1)) 
                       + I(log(hosted + 1)) 
                       + famine 
                       + I(log(battledeaths + 1)) 
                       + I(log(people_disaster + 1)) 
                       + trade_open 
                       + colony_top10 
                       + unsc_mtgs 
                       + any_pko 
                       + log(news) 
                       + as.factor (iso3) 
                       + as.factor(year), 
                       data = aid, 
                       family = "binomial" (link = "logit"))

summary (fit_logit_fe)

## END